/**
 * @file normal_distr.h
 *
 * The class to represent univariate normal distributions
 * 
 * @author Dahua Lin
 */

#ifndef SLIPP_NORMAL_DISTR_H_
#define SLIPP_NORMAL_DISTR_H_

#include <slipp/prob_base/pdistr_base.h>

namespace sli
{

	template<class RStream>
	class normal_sampler
	{
	public:
		SLIPP_DEFINE_SCALAR_DISTR_TYPES(double)

	public:
		normal_sampler(RStream& rstream, const scalar_type& mu, const scalar_type& sigma)
		: m_rng(rstream), m_mu(mu), m_sigma(sigma)
		{
		}

		sample_type draw()
		{
			return m_rng.get() * m_sigma + m_mu;
		}

		sample_collection_type draw(size_t n)
		{
			sample_collection_type v((index_t)n);
			m_rng.get_m(n, v.data());
			return v * m_sigma + m_mu;
		}

	private:
		normal_rng<RStream> m_rng;
		double m_mu;
		double m_sigma;

	};


	class normal_distr
	{
	public:
		SLIPP_DEFINE_SCALAR_DISTR_TYPES(double)

		class param_type
		{
		public:
			param_type() : m_mu(0), m_sigma(1) { }

			param_type(const scalar_type& mu, const scalar_type& sigma)
			: m_mu(mu), m_sigma(sigma) { }

			scalar_type mu() const
			{
				return m_mu;
			}

			scalar_type sigma() const
			{
				return m_sigma;
			}

			bool operator == (const param_type& rhs) const
			{
				return m_mu == rhs.m_mu && m_sigma == rhs.m_sigma;
			}

			bool operator != (const param_type& rhs) const
			{
				return !(operator == (rhs));
			}

		private:
			scalar_type m_mu;
			scalar_type m_sigma;
		};

	public:
		// constructors

		normal_distr()
		: m_param()
		{
		}

		normal_distr(const scalar_type& mu, const scalar_type& sigma)
		: m_param(mu, sigma)
		{
		}

		normal_distr(const param_type& param)
		: m_param(param)
		{
		}

		// parameters

		const param_type& param() const
		{
			return m_param;
		}

		scalar_type mu() const
		{
			return m_param.mu();
		}

		scalar_type sigma() const
		{
			return m_param.sigma();
		}

		// statistics

		double mean() const
		{
			return mu();
		}

		double var() const
		{
			return math::sqr(sigma());
		}

		double entropy() const
		{
			return 0.5 * (math::log(SLI_M_PI_2X * var()) + 1.0);
		}

	public:
		class pdf_evaluator_type
		{
		public:
			typedef normal_distr::sample_type sample_type;
			typedef normal_distr::sample_collection_type sample_collection_type;

			pdf_evaluator_type(const param_type& param)
			{
				mu = param.mu();
				scalar_type s2 = math::sqr(param.sigma());
				J = math::rcp(s2);
				C = math::log(SLI_M_PI_2X * s2);
			}

			SLIPP_ENSURE_INLINE double logpdf(const scalar_type& x) const
			{
				return (-0.5) * (J * math::sqr(x - mu) + C);
			}

			SLIPP_ENSURE_INLINE double pdf(const scalar_type& x) const
			{
				return math::exp(logpdf(x));
			}

			darow_t pdf(const sample_collection_type& X) const
			{
				return map(_detail::make_pdf_eval_functor(*this), X);
			}

			darow_t logpdf(const sample_collection_type& X) const
			{
				return map(_detail::make_logpdf_eval_functor(*this), X);
			}

		private:
			scalar_type mu;
			scalar_type J; 	// = 1 / sigma^2
			scalar_type C;	// = log(2 * pi* sigma^2)
		};

		pdf_evaluator_type get_pdf_evaluator() const
		{
			return pdf_evaluator_type(m_param);
		}

	public:
		// sampling

		template<class RStream>
		normal_sampler<RStream> get_sampler(RStream& rstream) const
		{
			return normal_sampler<RStream>(rstream, mu(), sigma());
		}

	private:
		param_type m_param;
	};


	// traits

	SLIPP_DEFINE_BASIC_DISTR_TRAITS(normal_distr, normal_sampler)


}

#endif 
